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A COMPUTERIZED SYMBOLIC INTEGRATION TECHNIQUE FOR 
DEVELOPMENT OF TRIANGULAR AND QUADRILATERAL 
COMPOSITE SHALLOW-SHELL FINITE ELEMENTS 

C, M. Andersen* and Ahmed K. Noor** 

Langley Research Center 

SUMMARY 

Computerized symbolic integration is used in conjunction with group -theoretic tech- 
niques to obtain analytic expressions for the stiffness, geometric stiffness, consistent 
mass, and consistent load matrices of composite shallow shell structural elements. The 
elements are shear flexible and have variable curvature. A stiffness (displacement) 
formulation is used with the fundamental unknowns consisting of both the displacement and 
rotation components of the reference surface of the shell. Both triangular and quadrilat- 
eral elements are developed. The triangular elements have six and ten nodes. The quad- 
rilateral elements have four and eight nodes and can have internal degrees of freedom 
associated with displacement modes which vanish along the edges of the element (bubble 
modes). 

The stiffness, geometric stiffness, consistent mass, and consistent load coefficients 
are expressed as linear combinations of integrals (over the element domain) whose inte- 
grands are products of shape functions and their derivatives. The evaluation of the ele- 
mental matrices is thus divided into two separate problems - determination of the coef- 
ficients in the linear combination and evaluation of the integrals. The latter problem can 
be computationally the more time consuming and the present study focuses on simplifi- 
cations and means of reducing the effort involved in this task. 

The integrals are performed symbolically by using the symbolic-and-algebraic- 
manipulation language MACSYMA. The efficiency of using symbolic integration in the 
element development is demonstrated by comparing the number of floating-point arith- 
metic operations required in this approach with those required by a commonly used 
numerical quadrature technique. 


*Senior Research Associate in Mathematics and Computer Science, College of 
William and Mary, Williamsburg, Virginia. 

**Research Professor of Engineering, Joint Institute for Acoustics and Flight 
Sciences, George Washington University, Hampton, Virginia. 



INTRODUCTION 


Although considerable literature has been devoted to the analysis of isotropic shells 
using doubly curved finite elements, investigations of the finite -element analysis of lami- 
nated composite shells are rather limited in extent. The reliable prediction of the re- 
sponse of composite shells often requires the use of high-order shear -flexible elements 
with the consequent increase in the developmental effort of the elemental matrices (the 
stiffness, geometric stiffness, and mass matrices). 

In most of the published work, numerical integration was used for the evaluation of 
the element stiffness matrices. Although numerical integration is both simple and adapt- 
able to computer programing, it can become computationally very expensive, particularly 
for higher order elements. This drawback has been recognized and improvements have 
been suggested (refs. 1 and 2); however, the difficulty has not been overcome. This raises 
the question as to whether an alternative approach such as using anal 5 d;ic (closed form or 
symbolic) integration in the element development is computationally more efficient. The 
present study addresses this question. More specifically, the objective of this paper is to 
demonstrate the computational advantages resulting from the use of computerized sym- 
bolic integration in conjunction with group -theoretic techniques (or symmetry transforma- 
tions) in the development of triangular and quadrilateral shallow shell elements. 

The analytical formulation is based on a form of the shallow-shell theory modified 
such that the effects of shear deformation and rotary inertia are included. Indlcial nota- 
tion is used throughout the development, since it is particularly useful in identifying the 
symmetries. The integrals are performed S 3 nubolically by using the symbolic-and- 
algebralc-manipulation language MACSYMAl (refs. 3 and 4). This is done in order to re- 
duce the tedium of analysis and to decrease the likelihood of errors. The use of group - 
theoretic techniques (or symmetry transformations) results in considerable reduction in 
both the amount of analytic computation needed and the size of the programs for numeri- 
cal computation. 


SYMBOLS 


basic integrals defined in equations (14) to (16) 
a*^,d*^,e*^ representations of dihedral groups defined in equations (56), (60), (61), and (69) 


^The MACSYMA system is being developed by the Mathlab group at Massachusetts 
Institute of Technology under the support of ARPA (Advanced Research Projects Agency 
of the U.S. Department of Defense) through Office of Naval Research Contract 
No. N00014-70-A -0362-001. 


2 



I 


[C] 8 by 8 matrix of shell stiffnesses 

extensional and transverse shear stiffnesses of shell, respectively 


^a^yp’^aZ^Z extensional and transverse shear stiffnesses of kth layer of shell 
C (7 portion of shell boundary over which tractions are prescribed 


^a^yp bending stiffnesses of shell 


EL,Erp 


®Q!/3 


elastic moduli in direction of fibers and normal to fibers, respectively 
permutation symbol 


' a fiyp 


stiffness interaction coefficients of shell 


cr^l 

'^a/30> '^a/31’ 

cr^j crij 
^a02>'^a/3Z 


functions of coordinates of corner nodes 


GlTjGtT shear moduli in plane of fibers and normal to plane of fibers, respectively 


Qa^O’Qor/Sl 

hlohk-1 

[K] 


functions of coordinates of corner nodes of elements 
weighting coefficients for numerical quadrature 

distances from reference (middle) surface to top and bottom surfaces of kth 
layer, respectively 

element stiffness matrix 

stiffness coefficients of shell element 
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[K] geometric stiffness matrix 

KjIj- geometric stiffness coefficients of shell element 

curvatures and twist of shell reference surface 

nodal values of 

L]^,L 2 ,Lj,L 2 ,L logarithmic functions defined in equations (33), (38), and (42) 

[M] consistent mass matrix 

mJj consistent mass coefficients of shell element 

^a/3 prescribed bending stress resultants on shell boundary 

mo,mi,m 2 density parameters of shell 

shape or interpolation function 

uAJq,^ prescribed extensional (in -plane) stress resultants on shell boundary 

relative magnitudes of initial stress resultants (prestress components) 
n^; unit outward normal to shell boundary 

consistent load vector 
Pj consistent load coefficients 

Pq,,P external load intensities in coordinate directions 

pjj,,p^ nodal values of p^, and p 
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2/3 


prescribed transverse shear stress resultants on shell boundary 


Rj coefficients associated with A-integrals 

integers associated with representative A-integrals 
r number of shape functions associated with an element 


ijk ijk ijk 

t>l ,b2 ,t>3 


coefficients associated with B- integrals 


integers associated with representative B-integrals 

s ratio, U2/U1 

s a quantity defined in equation ( 39 ) 

T kinetic energy of shell 

T transpose symbol 

coefficients associated with C-integrals 



integers associated with representative C-integrals 


t 

t 

U 

U° 


ratio, U3/U1 

a quantity defined in equation (39) 

strain energy of shell 

strain energy due to prestress 


Ul,U2,U3 


functions of nodal coordinates 


Uq!,W 


displacement components in coordinate directions (see fig. 1) 


= = = / linear functions of coordinates of corner nodes 

KvKzKzJ 

W work done by external forces 

X(j,X 3 orthogonal coordinate system (see fig, 1) 

Xq; values of x^, for corner nodes 

\ in-plane load parameter 

p number of numerical quadrature points 

^LT Poisson's ratio measuring strain in T direction due to uniaxial normal 

stress in the L direction 

local dimensionless coordinates in the shell domain 
I1 (uq,,w,<j!)q;) functional defined in equation (1) 
density of kth layer of shell 


a,T 

0Qf 

SI 


permutations 
rotation components 

nodal displacement vector 

nodal displacement parameters 

shell domain 

shell element domain 
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■■■■■I IBilBIHII II mill 


io circular frequency of vibration of shell 

= d/dXa 

The range of the different indices is as follows: 

Lower case Latin indices: 1 — r 

Upper case Latin indices: 1 — 5 

Greek indices: 1,2 

The finite -element models are designated: 

SQr stiffness formulation, quadrilateral element, r shape functions per funda- 

mental unknown 

STr stiffness formulation, triangular element, r shape functions per fundamental 

unknown 

The groups are designated: 

T )3 six-element dihedral group 

'T>4 eight -element dihedral group 

FINITE -EL EM ENT FORMULATION 

The analytical formulation is based on a form of the shallow-shell theory with the 
effects of shear deformation, anisotropic material behavior, rotary inertia, and bending- 
extensional coupling included, (See ref. 5.) The prebuckling state of the shell is assumed 
to be a momentless state. A displacement (stiffness) formulation is used in which the 
fundamental unknowns (dependent variables) consist of five generalized displacements — 
the two tangential displacement components uq., the normal displacement w, and the two 
rotation components cpa- (See fig, 1.) Throughout this paper the range of the Greek in- 
dices is 1,2 and it is implied that any index Greek or Latin, which appears twice in the 
same term is to be summed over. 
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(a) Displacements and loads. (b) Rotations. 

Figure 1.- Shell element and sign convention. 

Governing Functional of Shell 

The fundamental unknowns are assumed to vary sinusoidally with angular velocity 
w (the circular frequency of vibration of the shell). The functional used in the develop- 
ment of the elemental matrices is given (see ref. 6) by 

n(uQ(,w,0Q,) = U + U° - W - T (1) 


where U and U° are the strain energies 
potential energy of external forces, and T 

The expressions of U, U°, W, and 

are 


due to deformation and prestress, W is the 
is the kinetic energy of the shell. 

T in terms of the generalized displacements 


U = i ^ 8yUp + 2koip 9yUp w + kyp (w)^ + 2Fo;^ypj^8Q;U|3 8y^p 


+ kap dy4>p w + Da^yp By(f>p + 8q;W d0w + 2(f) a 9|3W + 4>a^^ | dO (2) 


U° = 



'^a/3 9q!W 8/3W dO 


W = \ (PofUcK + pw) 


4 . ,M. 


(3) 

(4) 
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and 


T = ^ y [mo(uQUo; + (w)2] + 2mi\ia<t>a + 


dS2 


(5) 


where ^a^yp> ^aj3yp> ^a^yp extensional, bending, and stiffness interac- 

tion coefficients of the shell; Ca2^2 are the transverse shear stiffnesses of the shell 
(see appendix A); are the curvature components and twist of the shell surface; 

are the initial stress resultants (prestress field), \ is the in-plane load param- 
eter; Pq, and p are the external load components in the coordinate directions Xq, and 
X3, respectively; mo, mj, and m2 are density parameters of the shell defined in 
appendix A; JU 0,^3, and are prescribed extensional, bending, and transverse 

shear stress resultants, respectively; is the shell domain; is the portion of the 
boundary over which tractions are prescribed; and no; is the unit outward normal to the 
boimdary. 


Finite -Element Discretization 

The shell domain is partitioned into triangular or quadrilateral subdomains 
With each is associated a finite element for which the fundamental unknowns 

are approximated by expressions of the form 


Uq, = 


w = NV3 ) 


( 6 ) 


where are the shape (or interpolation) functions; (J = 1 — 5) are displacement 

parameters; r equals the number of shape functions in the approximation; and a re- 
peated lower case Latin index denotes summation over the range 1 — r. Each set of dis- 
placement parameters is either associated with a node of the element or 

with a bubble function. 

The curvatures kQ,^ and the external load components Pq, and p are approxi- 
mated by the same shape functions used in approximating the generalized displacements, 
that is. 
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( 7 ) 


Pa = N^Pa > 
p = Nipi ^ 


where kj^| 3 , pj^, and p^ are the nodal values of p®, and p, respectively. To 

simplify the development of the elemental matrices, it is assumed in subsequent sections 
that the shell stiffnesses C’s, F’s, and D's as well as the density parameters mg, 
mj, m 2 , and the initial stress resultants are constants within each element. Equa- 
tions (2) to (5) can be esqjressed in terms of the shape functions and the nodal param- 
eters H/j as follows (for simplicity, the line integrals in eq. (r) are assumed to make no 
contribution) : 


u = i 


^ J (e) Ca,/3yp(9o!N^ + 2k”|3 


2 i-/ wov 

Elements “ 


k"^ k^p 8yNj + k^^ ayN^ 

^a/3yp 9aN^ 9yN^ *^+3 + ^a3/33 (9aN^ ^a+Z ’^S 


+ nIn3^J,^3 4^3) 


df2 


( 8 ) 


= 2 I £(e) 9 aNi P^xp{ df2 

Elements 


( 9 ) 


^ " Jj I (e) + P^4) 

ni j._ 


(10) 


Elements 


T = isl- 

2 ^ 
Elements 


1 Lie) N‘Nj[mo('('L'»4. + '''Ual 


dS2 (11) 
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The governing equations for each element are obtained by applying the stationary 
conditions for the functional in equation (1). This leads to the following set of equations: 




( 12 ) 


where kJJ are_ the stiffness coefficients; Kjj are the geometric stiffness (prestress) 
coefficients; Mjj and Pj are consistent mass and load coefficients; and X is an in- 
plane loading parameter used to determine the bucklii^ load or the magnitude of the pre- 
stress. The formulas for the aforementioned stiffness, mass and load coefficients are 
given in appendix B, For static stress analysis problems, \ = w = 0; and for free vibra- 
tion problems, \ = 0 and Pj = 0. The K*s, K’s, and M’s are completely symmetric 
under the interchange of one pair of indices for another, each pair of indices consisting of 
a superscript and the subscript just beneath it. 

To write equation (12) in matrix form, the first superscript -subscript pair of each 
of the K's, K's, and M's defines the row number and the second pair defines the col- 
umn number. For example, the term Kjj will be located in the [5(i - 1) + l]th row and 
in the [5(j - 1) + J]th column of the element stiffness matrix. Similarly, the P's and 
i/^'s can be written in column vector form. Thus, equation (12) becomes 


([K] + X[K]){^ = {P> + co2[m]{i/^> 


(13) 


Representation of Fundamental Unknowns 

Two triangular and four quadrilateral conforming elements were developed for the 
present study. The nodes of these elements are shown in figure 2, and their characteris- 
tics are summarized in table I. The triangular elements developed have six or ten nodes 
and the quadrilateral elements have four or eight nodes. For the two triangular elements 
and for two of the quadrilateral elements the range r of the superscripts in equations (6) 
and (7) is the number of nodes (that is, 6, 10, 4, or 8). For the two remaining quadrilat- 
eral elements, r is one greater than the number of nodes (that is, r = 5 or 9). The 
additional shape function, which vanishes along the entire element boundary and is not 
associated with a node, is usually referred to as a "bubble function." 

Since for each element all five of the fundamental unknowns are approximated by the 
same set of shape functions (see eq. (6)), the number of degrees of freedom per element is 
5r and thus for the elements developed the element stiffness matrices [K] of equation (13) 
range in size from 20 by 20 (for SQ4) to 50 by 50 (for ST 10). 
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(c) SQ8 and SQ9. (d) STIO. 

Figure 2.- The finite elements used in the present study. 


The nodal parameters introduced in equation (6) were selected to be values of 
the fundamental unknowns at the different nodes, and the approximations to the fundamental 
unknowns are continuous across the element boundaries. 

COMPUTATION OF ELEMENT CHARACTERISTICS 

Each component of the arrays Kjj, Kjj, Mjj, and Pj is a linear combination of 
integrals whose integrands are products of shape functions and their derivatives. (See 
appendix B.) If a constant shape function N® = 1 is introduced, then the integrals in 
appendix B can be reduced to linear combinations of integrals of the following three types: 


Aijkf ^ r NiN^N^N^ df2 

t 

O 

11 

(14) 

Bjj’^ = y (e) 

(i, j = 0 r; k = 1 - r) 

(15) 
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(i, j = 1 - r) 


(16) 


= y (e) dO 

The evaluation of the stiffness, mass, and load coefficients can be divided into two 
separate problems: (a) the determination of the coefficients in the linear combinations 
and (b) the evaluation of the A-, B-, and C -integrals. Since problem (a) is relatively 
simple, the present study focuses on simplifications and means of reducing the effort in- 
volved in problem (b). 


EVALUATION OF INTEGRALS 

For the evaluation of the A-, B-, and C -integrals (eqs. (14) to (16)), it is convenient 
to introduce a dimensionless local coordinate system (q! = 1,2) in terms of which the 

limits of integration are simply expressed. For triangular and quadrilateral elements, 
the 4 q! are chosen to be the area and quadrilateral coordinates, respectively. (See 
ref. 7 and fig. 3.) In the present study, the shape functions are simply polynomials in 
(see appendix C) and the Cartesian coordinates Xq, are linear or bilinear functions of 


\ 


I \ 




Figure 3.- Local coordinate systems for quadrilateral and triangular finite elements. 


The relations between the Cartesian and local dimensionless coordinates are 


dsi,£<!iM2) 

9 (^ 1 »^ 2 ) 


dll d|2 


(17) 
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and 


9 _ 9 

9 Xqj SXq/ B ^y/ 


9(xi,X2) 


9X|3 a 
®«/3 eyp ^ 


( 18 ) 


where 9(xi,X2y9(ll,l2) Jacobian determinant and d^yjdxQ^ is the inverse of the 

Jacobian matrix. See ref. 8 for more explicit expressions for 9(xi,X2)/9^|j,^2) 

9/9Xq; for selected cases. 

Through the use of equations (17) and (18), the A-, B-, and C -integrals can be ex- 
pressed in terms of the local coordinates |q, as 

Aijkf = y dll d|o (19) 


By>= = e„3e^py^(,)NlNi^|^d{,d£ 


( 20 ) 


- ^ay ^I3p 


r 9N^ 9NJ 
9|k 


9Xp 

9(xi,X2) 

9|j, 

9(1 1 , 12 ) 


-1 


dll d|2 


( 21 ) 


The A-integrals are completely symmetric under interchange of indices, and the 
B- and C -integrals satisfy the following symmetry relations: 




T^j i k 
= Ba 




j - 

^Q!/3 


- t^/SO! 


(22) 


Thus, the number of numerically distinct A-, B-, and C-integrals in equations (19) to (21) 
is (r+4)!/(4!r!), i(r+l)(r+2)(2r), and ■i(2r)(2r+l), respectively. (See ref. 9 and table n.) 

When the integrands of equations (19) to (21) are polynomials in |q,, the analytic 
(symbolic) evaluation of the corresponding integrals is a straightforward task. Examina- 
tion of equations (19) and (20) reveals that this is the case for all the A- and B-integrals. 
The integrands of the C-integrals (eq. (21)) are polynomials in |q, only when the Jaco- 
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bian determinant 9(xi,X2)y^9 independent of and this is the case for tri- 

angular and parallelogram elements. For a general quadrilateral element, the Jacobian 
determinant is linear in and consequently, the C -integrals are complicated functions 
of the nodal coordinates and involve logarithms. 

The forms of the analytic expressions for the A-, B-, and C -integrals for both the 
triangular and quadrilateral elements are given subsequently. 

Triangular Elements 

For triangular finite elements, the A-, and C -integrals may be expressed as 




(23) 





al 


V 


a2 


(24) 


Q 

II 


V 0 - 2 ] 


'T'ij 

^2 

^^1 




rpij 

/3 

rpij 

^4 

>2 


(25) 


where the R's, S's, and T's are rational numbers (independent of the nodal coordi- 
nates) and are discussed later; the Vs are functions of the nodal coordinates given by 


V 


al ~ 



Vq; 2 = e 




( 26 ) 


The quantity Uj is the area of the triangular element and is given by 


Ul = I Vq;iV|32 (27) 

12 

and x^, x^, and are the coordinates of the corner nodes numbered as in figure 2(b), 

Since the five quantities Uj, Vq,i, and Vq ,2 are functions of the nodal coordi- 
nates, they need to be evaluated for each element. On the other hand, the R's, S's, and 
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T's are the same for all triangular elements having the same number of nodes. As will 
be demonstrated in subsequent sections, the evaluation of the A-, B-, and C-integrals 
using eqiiations (23) to (25) is computationally more efficient than using a numerical quad- 
rature technique wherein the integrands have to be evaluated at each of the quadrature 
points, multiplied by the weighting coefficients, and then summed. 


Quadrilateral Elements 

For quadrilateral elements, the A-, B-, and C-integrals may be expressed as 






ijkf j^ijkf 




Ul 

U2 

U3 


(28) 


B 


ijk_ 


a 


4jk qijk _ijk 

>1 S2 S3 I 


V«l 


'of2 


V, 


a3 


(29) 


^a(3 - \^Q'/30 + 




9^ 


1] 

Qf/32 




ij 

q;/S3 


Li(s,t) 

L2(s,t) 

L2(t,s) 


(30) 


where the R’s and S's are rational numbers (independent of the nodal coordinates) and 
are discussed later; the U's and Vs are functions of the nodal coordinates given by 


Vai = 

Va2 = 
Vq-3 = 


^ais{ 4) 

4 

4 

^a^{-4 - 4 ^4-^4) 

4 


(31) 
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and 


Ul = Vq,2V|33^ 


U2 = ea^ Vq,3Vj3i V 


( 32 ) 


U3 = eo!|8 Vq,iV^2J 


1 9 o A 

and x^, x^, x^, and x^, are the coordinates of the corner nodes numbered as in fig- 
ure 2(a). The quantity Uj is equal to one -fourth the area of the element. 

The ‘3T's in equation (30) are complicated rational functions of s, t, and the Vs, 
and the L's are logarithmic functions given by 


Lj(s,t) = L (t,s) = I log 


1 - (s + t)^ 


L 2 (s,t) = I log 


1 - (s - t)2 

(1 + t)^ - 


(1 - t)2 - s2 




where s and t are dimensionless parameters given by 







(33) 


(34) 


The relations between the Cartesian coordinates Xq, and the local coordinates in- 

volve the U’s and V*s of equations (31) and (32). They are 


- Vq;1^1^2 


+ Vq;2^1 + ^aZ^2 


x^ + x^ + x^ + x^ 
4 


(35) 
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and 


= Ui - U2«2 - U3«l 


(36) 


The numerical evaluation of the A-integrals and B -integrals by using equations (28) 
and (29) is a straightforward task. However, since the Q^’s in equation (30) typically 
become infinite in the limit as s and/or t goes to zero, and since elements with small 
or vanishing s or t must often be evaluated, equation (30) is not a satisfactory form to 
use for the numerical evaluation of the C -integrals. 


For the C -integrals used in this study, the difficulties are overcome by casting 
equation (30) in the form 


■'Q'/S 


Ul 




Li(s,t) 


‘^a^l '^a^2 '3'’ffi33 

L2(s,t) \ 



L2(t,s)_ 



J 


(37) 


where 


Li(s,t) = Li(t,s) = (it) 




Li(s,t) + 2st + 



L2(s,t) = (t)-^ 


L2(s,t) - 2t - (I + 2s2jt 3 


(38) 


and 


s = 


1 - t" 


t 


t 

1 -s2 


(39) 


Since only convex quadrilaterals (those for which the Jacobian (36) does not vanish 
anywhere within the domain of the elements) are of interest, equation (37) needs to be 
evaluated only for r and s satisfying the condition 

|s| + |t| < 1 (40) 
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(Elements with |s| + |t| = 1 have three of their four vertexes lying on a straight line; 
that is, they have degenerated into triangles.) The ?T's are free of singularities except 
at s = ±1 and t = ±1, and their numerical evaluation is straightforward. Elements with 
s or t equal or nearly equal to one are not of interest anyway. 

Even though the functions Lj and L 2 are free of singularities within the region 
specified by equation (40), they must be evaluated with care. Equation (38) is used di- 
rectly for numerical evaluation only if neither s nor t is small; otherwise, Taylor 
series expansions are used. These series expansions were derived by using MACSYMA 
and are part of the computer program listed in reference 11. 

Although the A-, B-, and C -integrals for any quadrilateral element can be numeri- 
cally evaluated by using equations (28), (29), and (37), simplified forms for these integrals 
should be used to improve computational efficiency whenever the element is a trapezoid or 
a parallelogram. 


Trapezoidal Elements 

Trapezoidal elements are characterized by the vanishing of t = Us/Uj or else 
s = U 2 /U 1 . In the case t = 0, equation (37) reduces to 


C 


ij 

a/3 


_ 1 _ 

Ul 


Qa/SO + 9a/31 L(s) 


(41) 


where the 


9's are functions of s and the V’s, and 


L(s) 



-2s - f s3 

O 


(42) 


Similar expressions hold for the case s = 0. A Taylor series expansion is used to eval- 
uate the function L when its argument is small. 


Parallelogram Elements 

Parellelogram elements (and thus rectangular elements) are characterized by the 
vanishing of both s and t. The formulas for the A-, B-, and C-integrals over quadri- 
lateral elements (eqs. (28) to (30)) reduce in the parallelogram case to 
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(45) 


where the R's, S's, and T’s are rational numbers. These expressions, like those for 
the triangular elements (eqs, (23) to (25)), are entirely free of logarithms. 

GROUP-THEORETIC TECHNIQUES 


The evaluation of the A-, B-, and C -integrals can be divided into two distinct phases: 
(1) obtaining analytic (or algebraic) expressions for the integrals using the symbolic-and- 
algebraic-manipulation language MACSYMA and (2) numerically evaluating the resulting 
expressions. These two phases will henceforth be referred to as the symbolic and numer- 
ical phases, respectively. 

if k 

Because of the way in which the integrals B<;x depend on the values of the Greek 
index a (see eq. (20)), separate symbolic integrations for fixed i,j,k are not required 
for Bj^'^ and Similarly, for given i,j only one symbolic integration is required 

for the set of integrals C^^, C ^21’ ^22* consideration and by 

taking into account the symmetries of the Latin indices in equations (19), (20), and (21), 
the number of symbolic integrations required to determine the coefficients in the expres- 
sions (eqs. (23) to (25) and eqs. (28), (29), (37), (41), and (45)) for the A-, B-, and C- 
integrals would be (r+4)!/(4!r!), -|(r+l)(r+2)(r), and ^(r)(r+l), respectively. (The Q’s 
in eq. (41) and the T's in eq. (45) do not require separate symbolic integrations since 
they can be determined from the ^'s of eq. (37)). Thus, for the six elements developed 
in this study, a total of 4647 symbolic integrations would be required. The evaluation of 
such a large number of integrals symbolically, some of which are very complicated, seems 
discouraging even with the aid of a computer, and it is highly desirable to reduce this num- 
ber. This reduction can be accomplished for the triangular and quadrilateral elements 
studied herein by taking advantage of the similarities in algebraic form among the various 
integrals, viz, by employing the theory of group transformations (refs, 9 and 10). 
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Basic Idea of Group -Theoretic Techniques 

The essential idea of the group -theoretic techniques employed herein is that a rela- 
tively small number of the A-, B-, and C -integrals are selected as "representative^' inte- 
grals for which symbolic calculations are performed directly. Each of the other A-, B-, 
and C -integrals is related to one of these representative integrals by symmetry consider- 
ations. Therefore, the symbolic phase of computation reduces to obtaining analytic ex- 
pressions for just the representative integrals. The numbers of representative integrals 
required for each of the various elements considered in the present study is given in 
table II. The total number of representative integrals is 886; thus, the group -theoretic 
techniques reduce the number of symbolic calculations to less than one -fifth of what it 
would otherwise have been. 

In this study, group transformations are applied by using two different (but mathe- 
matically equivalent) procedures.* The first procedure is used for the sets of integrals 
that can be expressed as linear combinations (with numerical coefficients) of a small num- 
ber of functions of the nodal coordinates. In this case the output from the symbolic phase 
of computation consists of arrays of integers obtained from the symbolic integration of 
representative integrals. In the numerical phase the coefficients for all the integrals re- 
lated to these representative integrals are numerically evaluated and stored for repeated 
subsequent use. As a result, the evaluation of these integrals is very fast and only a rel- 
atively small amount of data needs to be transferred from the symbolic phase of computa- 
tion to the numerical phase. This procedure was found to be the more efficient one for 
evaluating the A- and B -integrals (for all elements studied herein) and the C -integrals for 
the triangular and parallelogram elements. However, it would be very cumbersome to try 
to apply this procedure to the remaining C -integrals. 

For evaluating the C -integrals for nonparallelogram quadrilateral elements, the sec- 
ond procedure involving group transformations was used. Since these integrals have more 
complicated forms, the output from the symbolic phase of computation consisted of alge- 
braic expressions for the representative integrals rather than just arrays of coefficients. 
These algebraic expressions are used for the numerical evaluation of the nonrepresenta- 
tive as well as of the representative C -integrals. This is done by applying suitable group 
transformations to the variables contained in these expressions before evaluating them. 

By using the same algebraic expressions for several related integrals, fewer symbolic 
calculations are performed and computer memory requirements for the numerical phase 
of computation are reduced. 

The second procedure could have been applied to all the integrals but with some loss 
in the computational efficiency. The difference in the way the two procedures would apply 
to the A- and B-integrals is that by using the first procedure, the coefficients (the R’s 
and S's) are transformed (and this needs to be done only once), and by using the second 
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procedure the functions of the nodal coordinates (the U's and V*s) would be transformed 
(and this would have to be done for each element). 

Before discussing the application of these two procedures in detail, the concept of 
dihedral groups is briefly reviewed. 


Dihedral Groups 

For triangular and quadrilateral elements the symmetry groups which serve to re- 
late the various A-, B-, and C -integrals to their representative integrals belong to an in- 
finite family of groups called dihedral groups. (See ref. 10.) For triangular elements 
the relevant group is the six-element dihedral group (denoted by ^ 3 ), and for quadrila- 
teral elements it is the eight -element dihedral group (denoted by ^ 4 ). The group ^>3 
is the symmetry group of an equilateral triangle, and the group ^>4 is the symmetry 
group of a square. These groups are relevant because any triangular element can appear 
as an equilateral triangle when viewed in an appropriate local coordinate system, and any 
quadrilateral element appears as a square as viewed in the quadrilateral coordinate sys- 
tem of figure 3(a). 

To each element in the symmetry group there corresponds a linear transformation 
on the local coordinate system which maps the element domain onto itself. Let such a 
transformation be symbolized by 

= ( 46 ) 

This transformation induces a transformation on any function f(^i,^2) sense that 

(4'7) 


where 


In particular, the set of shape functions i = 1 "" r, has the property that 

each transformed shape function (n^) is identical to one of the untransformed functions 
Nl. In fact 


(n^)' = N^(‘) 


(49) 
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where r is a permutation on the set of r indices. Since the indices on the shape func- 
tions, in effect, form the superscripts on the A’s, B's, and C’s, the effect of the group 
transformations is to induce transformations on these superscripts. 

To each element in the symmetry group, there also corresponds a permutation a 
on the coordinates of the corner nodes 


X 


i 

a 




for triangular elements;^ 
quadrilateral elements j 


(50) 


These permutations induce linear transformations on various sets of functions of the x^. 
Examples of such sets are U^; U 2 and U 3 ; Vq,j^; Vq ,2 and Vq, 3 ; r and s; Lj; 
L 2 and L 3 ; Li; L 2 and L 3 ; etc. Thus, for example, 


Ui(xJk) Ui(xg(i)) = aUi(xiQ,) 

and 






U2(4) 




= d 


«s(4) 



1 

Us (4) 


(51) 


(52) 


where a is +1 or -1 and d is a 2 by 2 matrix. Both a and d depend on the particu- 
lar group transformation. In this way certain transformations on the superscripts of the 
aU«, Bjil', and are related to linear transformations on the variables contained 

in the formulas for these integrals. 

First Procedure Based on Group Transformations 

As mentioned previously, group transformations are employed in two different ways. 
The first procedure based on group transformations is used to evaluate the R’s, S’s, 
and T’s in equations (23) to (25), (28) to (29), and (43) to (45). 

For triangular elements the formulas used to evaluate these coefficients are 

( 53 ) 
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and 
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T§ 



( 54 ) 


(55) 


where the d’s are matrices which represent group transformations. The d's are 
given by 



(56) 


In equations (53) to (55) T denotes transpose, and the ^'s, S's, and ‘ST's are 
integers obtained from the symbolic integration of representative integrals. 

For quadrilateral elements the formulas used to evaluate the R’s, S’s, and T's 

are 


Rf' = (■'Rj')'’*!’ 
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where the a's and d's represent group transformations and are given by 


( 59 ) 


al = -a^ = a^ = -a^ = a® = -a® = aJ = -a^ = 1 (60) 
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d^ = -d^ 

d4 = -d2 
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'l 

0 
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-1 

d6 = 

’ 0 
-1 

0 

d7 = -d5 

d^ = -d® 


and the '15's, 2's, and ^'s are again integers obtained from the symbolic integration 
of representative integrals. 

In equations (53) to (55) and (57) to (59) the superscripts m, in, m, n, n, and h 
in the right-hand members must be understood as integer -valued functions of the super- 
scripts in the corresponding left-hand members, that is 


m = m(i,j,k,e) 


m = m(i,j,k) 


m 


= m(i,j) 


n = n(i,j,k,0 


n = n(i,j,k) 


n = n(i,j) 


(62) 


The superscripts m, m, and in indicate which representative integral is being refer- 
enced, and the superscripts n, h, and h indicate which group transformation is being 
applied. The dependence of the superscripts m, in, in, n, n, and h on the super- 
scripts i, j, k, and £ is given by a set of integer -valued input arrays generated by 
separate computer programs. These arrays are given in reference 11. 

Whenever the left-hand member of one of the equations (53) to (55) or (57) to (59) is 
one of the representative integrals, then the superscript n or n or H (as the case 
may be) equals one which indicates the identity group transformation. Thus for triangular 
elements the coefficients for the representative A-, B-, and C-integrals are 
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S5^/s“ and and 0 ’“/' 0 r 5 ^, 9'3^/'3'5^, and 'ET^/gr™; respectively; 

and similarly for quadrilateral elements. The arrays of <12-, S-, and ET -values are given 
in reference 11. 


First Procedure Illustrated 

To illustrate the use of this procedure, consider the integral for the ST4 

element. Since was selected as one ofthe representative integrals, it has n = 1. 

From the way in which the representative integrals have been ordered, A^^^^ is the 
eighth representative A-integral, and thus it has m = 8. Because A^m is a represen- 
tative integral, it is one of the integrals to be symbolically integrated. The result of this 
integration is 


A 


2111 = -Ui - 2U2 + 3U3 
75 


(63) 


and thus, by comparing equations (57) and (61) with equation (27), it follows that 


<^1 = -1 = - 2 ^ 

> 

<12| = 3 ^4 = 75 


(64) 


These four values are included in the <12 -array for the ST4 element assembled during 
the symbolic phase of the computation. 

Seven other A-integrals for the ST4 element are also characterized by m = 8. 
They are a 2221^ ^3222^ ^3332^ ^4333^ ^4443^ and A^m, and their n 

values are 5, 4, 6, 2, 8, 3, and 7, respectively. Thus, from equations (52) and (56), it 
follows that their coefficients are 




r 2111 = _o.0133 

r|111 = -0.0267 

= 0.0400 

R^Zl _o.oi 33 

r2221 ^ .0.0267 

r2221 ^ .0.0400 

Rp22 ^ _o 0133 

r|222 ^ .0.0400 

Rp22 ^ .0.0267 

r^ 322 ^ _o 0133 

r|^^2 ^ .0.0400 

= 0.0267 


j 


(Equations continued on next page) 
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r|333 ^ .0,0133 

r|^^^ = 0.0400 

r4333 ^ 0.0267 

r|443 ^ .0.0133 

rI***^^ = 0.0400 

r|443 ^ .0.0267 

r4441 ^ .0.0133 

= 0.0267 

r|441 ^ .0.0400 1 

r|111 = .0,0133 

r|^^^ = 0.0267 

r|^^^ = 0.0400 ^ 


These results are included in the R -array which is evaluated in the numerical phase of 
computation and stored for subsequent use. Then for each finite element processed, the 
U's are evaluated by using equation (30) and after this evaluation the set of distinct A- 
integrals (that is, the with i^j Sk? £) is numerically evaluated and saved in an 

array. Included in this A-array are the eight integrals A^Hl, A^221^ ^3222^ ^3332^ 
a 4333^ A'^443^ ^4441^ and a 41H, whose coefficients are given in equations (65). Fin- 
ally, as various A-integrals are needed in the evaluation of the element characteristic 
matrices, they are quickly retrieved from the A-array. Because of the permutational 
symmetry of the indices, the 32 A-integrals on the following list all depend on the ^- 
values given in equations (64) : 


^2111 ^ ^1211 ^ ^1121 ^ ^1112 
^3222 ^ ^2322 ^2232 ^ ^2223 

^4333 ^ ^3433 ^ ^3343 ^ ^3334 
^1444 ^ ^4144 ^ ^4414 ^ ^4441 


^1222 ^ ^2122 ^ ^2212 ^ ^2221 
^2333 ^ ^3233 ^ ^3323 ^ ^3332 
^3444 ^ ^4333 ^ ^4434 ^ ^4443 
^4111 ^ ^1411 ^ ^1141 ^ ^1114 


In this way there can be as many as 192 A-integrals associated with a single representa- 
tive A -integral. 


Second Procedure Based on Group Transformations 

The second procedure based on group transformations is used for the evaluation of 
C -integrals for nonparallelogram quadrilateral elements. Because of the greater com- 
plexity of the expressions for these integrals, the output from the symbolic phase of com- 
putation consists of FORTRAN-coded formulas for the representative integrals rather 
than arrays of integers as in the other cases. These formulas are functions of the vari- 
ables Ui, s, t, L]^(s,t), L 2 (s,t), L 2 (t,s), and the Vs. The procedure by which both 
a representative C -integral and the integrals related to it by dihedral symmetry can be 
evaluated by use of many of the same FORTRAN statements is the following; 
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(1) Based on the superscripts I and j, the program looks up the values of 

m = m(i,j) and n = n(i,j), where m designates the representative integrals and h des- 
ignates the group transformation. 

(2) The quantities Vq,i, Vq,2> Vq,3, s, t, L 2 (s,t), and L 2 (t,s) undergo the fol- 
lowing transformations: 


V«l ■* 


V«2 


V^2 


V«2 




= a" d" 


_V«3 

i 

vis 


>3 



( 66 ) 


(67) 


L2(s,t)' 


L§(s,t) 


L2(s>t) 

L2(t,s) 


L|(t,s) 

= 

L 2 (t,s) 


( 68 ) 


where 



(69) 


and the a's and d's are given in equations (60) and (61). (Computationally, it is much 
faster to use eqs. (65) than to reevaluate the logarithmic functions each time s and t 
are transformed.) The variables Uj and L|(s,t) = Lj(sa,ta) are not transformed. 
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(3) The FORTRAN code for the mth representative integral is executed by using 
these transformed variables in place of the original variables. This technique of using 
the same FORTRAN statements to evaluate several symmetry -related integrals signifi- 
cantly reduces the program length, 

PERFORMANCE EVALUATION AND COMPARISON 
WITH NUMERICAL QUADRATURE 

In order to assess the computational efficiency of the symbolic integration approach 
for developing elemental matrices, a computer program called SYMINSE (SYMb olically 
INtegrated Shell dements) was developed based on this approach. The program generates 
the element characteristic arrays (stiffness, geometric stiffness, mass, and load matrices) 
for the six elements in table I. A detailed description of SYMINSE program is given in 
reference 11. 

The CPU times required for the SYMINSE program to generate the element charac- 
teristic arrays for 14 different shallow-shell elements and 14 flat-plate elements (with 
bending extensional coupling included) are given in table III. Also included in this table 
are the field lengths required for these elements. No comparison is given of the CPU 
times required for computing the characteristic arrays for these elements via the symbolic 
integration approach compared with the numerical quadrature approach. This is because 
any such comparison would, by necessity, be influenced by the efficiencies of both pro- 
grams and would be dependent on compilers, operating systems, and computer hardware 
as well. 

Instead a quantitative measure of the efficiency of the symbolic integration approach 
is provided by giving a comparison of the number of floating-point multiplications and 
additions in a conventional numerical quadrature approach with the number required by 
the SYMINSE program. 

In the conventional numerical quadrature approach, the element stiffness matrix is 
expressed in the form 

£=1 

where [Bg] is the strain -displacement matrix for interpolation models, [C] is the matrix 
of shell stiffnesses, are the weighting coefficients for numerical quadrature, ju is 

the number of quadrature points, and r is the number of shape functions. 
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To generate the element stiffness matrix [K], the matrix of shell stiffnesses [C] is 
multiplied in turn by each of the weight coefficients Hjg. Because of the symmetry of the 
element stiffness matrix, only the upper triangvilar part of [K] needs to be evaluated. The 
individual contributions at each of the quadrature points are added to give the stiffness 
matrix. 

Table IV gives the nximber of primary floating-point arithmetic operations required 
to generate the matrix [K] as a function of the number of numerical quadrature points fx 
and the number of shape functions r. The numbers given in this table are based on the 
assumption that each element of [B{] is a single quantity (rather than an expression in the 
local coordinates) and therefore, the operation counts given are conservative, that is, less 
than the actvial number of arithmetic operations. 

The number of floating-point arithmetic operations required to generate the element 
characteristic arrays for triangular and parallelogram elements using the SYMINSE pro- 
gram are given in table V. Note that very few additional operations are required for the 
evaluation of the geometric stiffness, the consistent mass, and the consistent load coeffi- 
cients. The trapezoidal and nontrapezoidal quadrilateral (trapezium) elements are not in- 
cluded in table V because their arithmetic operation counts are more difficult to obtain. 
However, the CPU times listed in table III give a measure of the additional complexity for 
these cases. 

As a test for the reliability of the operation count estimates, the ratios of the CPU 
times expended in evaluating integrals compared with the CPU times expended in forming 
the characteristic arrays as linear combinations of integrals are given in table VI. These 
ratios may be compared with the comparable ratios of operation counts also given in 
table VI. The ratios are not very close but agree in order of magnitude. 

Figure 4 shows the ratios of the arithmetic operation counts involved in the evalua- 
tion of the element stiffness array for both the symbolic manipulation and numerical quad- 
rature procedures. 

The element characteristic arrays computed by the SYMINSE program have been 
verified by comparing them with those obtained by a program based on numerical quadra- 
ture. This latter program, although not optimized for efficiency, ran much slower than 
the SYMINSE program. This fact together with the preceding comparison of operation 
counts strongly suggest that the symbolic integration approach presented herein holds 
much promise. 


POSSIBLE FURTHER IMPROVEMENTS AND EXTENSIONS 

On the basis of the present study, the following remarks regarding the use of sym- 
bolic integration in the development of higher order elements seem to be in order: 
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Number of quadrature points 



Number of quadrature points 

Figure 4.- Comparison of number of floating-point arithmetic operations required 
to evaluate integrals using numerical quadrature and symbolic integration tech 
niques. r, number of shape functions. 






(1) The analytic methods discussed herein can be easily implemented without modi- 
fication for triangular elements with more than ten nodes and for parallelogram elements 
with more than eight nodes, 

(2) A minor modification of these methods is required for trapezoidal and trapezium 
(nontrapezoid quadrilateral) elements with more than eight nodes. For these elements 
the 'Sf's and Q's of equations (37) and (41) contain singularities when s or t equals 
zero. These difficulties are eliminated by using logarithmic functions obtained by sub- 
tracting off low order terms in the Taylor series expansions of L(s), Li(s,t), and 
L 2 (s,t) defined in equations (42) and (38). 

(3) A minor modification of the way in which symmetry transformations are applied 
is required for quadrilateral elements with more than one bubble mode. This case is dis- 
cussed in reference 9. 

(4) With these modifications the implementation of trapezoidal elements with more 
than eight nodes and/or with more than one bubble mode is not difficult. However, the 
implementation of trapezium elements with more than eight nodes or with more than one 
bubble mode would be difficult because of the complexity of the analytic expressions, 

(5) For the trapezium elements a hybrid approach which combines the exact symbolic 
integration with either approximate anal 5 dic integration or numerical quadrature appears 
to have advantages over the purely numerical quadrature and the purely symbolic integra- 
tion approaches. In the hybrid approach the A- and B-integrals are evaluated exactly by 
using symbolic integration and the C -integrals are evaluated approximately by numerical 
quadrature or symbolic integration. This approach would retain the major advantages re- 
sulting from symbolic integration of the A- and B-integrals and eliminate the difficulties 
associated with exact symbolic integration of the C -integrals. 

A count of floating-point arithmetic operations suggests that even for a purely 
numerical quadrature approach, the evaluation of A-, B-, and C-integrals followed by 
formation of the stiffness coefficients as a linear combination of these integrals is faster 
than the conventional approach discussed in the preceding section. This suggests that 
symbolic integration and numerical quadrature can be readily combined in one program, 

CONCLUDING REMARKS 

Computerized symbolic integration and group -theoretic techniques are employed in 
combination to obtain analytic expressions for the stiffness, geometric stiffness, consis- 
tent mass and consistent load coefficients of composite shallow-shell structural elements. 
The elements are shear flexible and have variable curvature. A stiffness (displacement) 
formulation is used with the fimdamental unknowns consisting of both the displacement and 
rotation components of the reference surface of the shell. Both triangular and quadrilat- 
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eral elements are developed. The triangular elements developed have six and ten nodes. 
The quadrilateral elements have four and eight nodes and can have internal degrees of 
freedom associated with displacement modes which vanish along the edges of the element. 

The coefficients of the elemental matrices (stiffness, geometric stiffness, consistent 
mass, and consistent load) are expressed as linear combinations of integrals (over the 
element domain) whose integrands are products of shape functions and their derivatives. 
The evaluation of the elemental matrices is thus divided into two separate problems — 
determination of the coefficients in the linear combination and evaluation of the integrals. 
The latter problem can be computationally the more time consuming and the present study 
aims at reducing the computational effort involved in this task through the use of compu- 
terized symbolic integration. 

Based upon a comparison between the number of floating-point arithmetic operations 
required in the symbolic integration and traditional numerical quadrature approaches, it 
appears that the symbolic manipulation approach, as presented herein, is considerably 
more efficient and holds much promise. This is especially true for the triangular, paral- 
lelogram, and trapezoidal elements. 

For trapezium (nontrapezoidal quadrilateral) elements, however, a hybrid approach 
in which some of the integrals (the A- and B -integrals) are evaluated exactly by symbolic 
manipulation and other integrals (the C -integrals) by either numerical quadrature or by 
means of approximate analjd;ic expressions appears to have advantages over a purely 
numerical quadrature or a purely symbolic integration approach. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, Va. 2366 5 
October 3, 1975 
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APPENDIX A 


ELASTIC COEFFICIENTS OF LAMINATED SHELLS 


Elastic Stiffnesses of Layers 

The nonzero stiffness coefficients and of the kth orthotropic layer 

of the shell referred to the directions of principal elasticity are given by 


^1111 


„(k) 

^^1122 


(k) 

^=2222 


c(k) 

*^1212 


c(l^) 


^2323 







where the subscripts L and T denote the direction of fibers and the transverse direc- 
tion, respectively, the superscript k refers to the kth layer, the E’s are elastic mod- 
uli, the G's are shear moduli, and is the Poisson’s ratio measuring the strain in 

the T -direction due to a uniaxial normal stress in the L -direction 


'^TL^L = i'LT^'T 


\ = 1 _ i^lt^TL 


The stiffness coefficients Co-/3yp and satisfy the following symmetry 

relationships: 

^a^yp ~ ^ypa^ ~ ^^ayp ~ ^afipy 
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and 


Cq!3^3 = C|33q!3 = C3q!| 33 = Ca33/3 

If the coordinates Xq, are rotated, the elastic coefficients CQ,^yp and 
transform as components of fourth- and second-order tensors, respectively. The trans- 
formation law of these coefficients is expressed as follows: 

^a'jS'y'p' “ ^oi^YP ^oioi' Vx' ^P' 


and 


Cq!'3^'3 “ Cq;3^3 laa^ 

where and Cq,' 3 | 3'3 are the elastic coefficients referred to the new coordi- 

nate system Xq,, and 

iaa' = cos (xq;,Xq,') 


Elastic Coefficients of Shell 

The equivalent elastic stiffnesses of the shell are given by 


NL 


'afiyp ^a(3yp ^al3yp, 


"hk (k) 
^a(3yp 

k=l ^k-1 




1 Xc 


x| dx «3 


and 


NL 


„ V (k) , 

Cq-3/33 = ^ J Cq, 3^3 dxg 

k=l ^k-1 


where NL is the total number of layers of the shell and and are the dis- 

tances from the reference surface to the top and bottom surfaces of the kth layer, respec- 
tively. The elastic compliances of the shell ^alBypy ^ajBypy ^a3l33 

obtained by inversion of the matrix of the elastic stiffnesses. (See ref, 5.) 
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The shell stiffnesses and compliance coefficients satisfy symmetry and transforma- 
tion relations similar to those of the stiffness coefficients of individual layers. 

The density parameters of the shell are given by 



where is the mass density of the kth layer of the shell, 

O 


36 


j 



■ ■■!■ II Mil III I III I II I 


II III 


APPENDIX B 

FORMULAS FOR THE COEFFICIENTS IN GOVERNING EQUATIONS 
FOR INDIVIDUAL ELEMENTS 

The expressions for the independent stiffness coefficients in equation (12) are given 

by 

Kq!( 8 = ^oCe) 

' L(e) n"n' an 

Ka,j+3 = f (e) ^ a^rp ^P^^ 

Kg^ = y (g) {CaY^P kay kg, N%jN%”^ + C 0,3^3 OaN^ d^N^) 

4,L+3 = y (e) (Fo-y/Sp 4 p ^^N^ + 0^3/33 N^ a^N^) dSl 

Ka+3,0+3 = y (e) (Day/Sp ^yN^ 9pN^ + dS2 

The independent nonzero geometric stiffness coefficients are given by 
^33 = y (e) 

a 


The independent nonzero consistent mass coefficients are given by 

Kh^U) ^0 N%j dO 

06' 

M^i+3 = y^(e) “1 dS2 
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moNW dJ2 

Ma+3 J+3 = f (e) “^2 ^a/3 dfi 
0^' 


where Sq,^ is the Kronecker delta on a and 

Finally, the consistent nodal load coefficients are given by 

N'N^pm dS2 
’ In(e) n'n^p"” dn 

In these equations the range of the lower case Latin indices is 1 — r where r is 
the number of shape functions, the range of the Greek indices is 1,2, and a repeated index 
denotes summation over the full range of the index. However, for the elements with inter- 
nal degrees of freedom (SQ5 and SQ9) the indices S. and m in these expressions have 
a range equal to the number of nodes in the element (that is, 4 for the SQ5 element and 
8 for the SQ9 element). This amounts to distributing the loading on the nodes of these 
elements with no loading associated with the internal degrees of freedom. 
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APPENDIX C 


SHAPE FUNCTIONS USED IN THE PRESENT STUDY 

The shape functions used In the present study are numbered as shown in figure 2. 

Four -Node Quadrilateral Elements (SQ4 and SQ5) 

The shape functions for SQ4 are 

4 

4 

4 

1.^4 - ^ ^2) 

4 

The shape functions for SQ5 are these functions plus the additional function representing 
the bubble mode. 

n5= (1 -4i2)(i -l2^) 


Eight -Node Quadrilateral Elements (SQ8 and SQ9) 
The shape functions for SQ8 are 

, (1 - {l)(l - S 2 X-I - ^1 - i2) 

4 

„i, (1 4 il)(l - l2)(-l ^ h - h) 

4 

(1 S2K-1 * li * ia) 
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. t 4 ^2)(-l - ^1 + ^2) 

^ “ 4 

2 

2 

(l-jl^)(l^t2) 

2 

„« (1 ^ {l)(l - i2^) 

2 

For SQ9 the additional shape function N® is 

n 9 = (1 - {i2)(l - {32) 

Triangular Elements (ST6 and ST 10) 
The shape functions for ST6 are 

n 1 = li(25i - 1) 
n2 = {2(252 - 1) 

n 3 = (1 - {j - {2)(1 - 2{i - 253) 

= 4{il2 
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n5 = 4{2(1-?i-?2) 

N« = 4 |i(l -^1 - ^2) 



aiiH III lull III 


f} 

I 

j 


APPENDIX C 

The shape functions for STIO are 

^l(34i - l)(3|i - 2) 

2 

^2 _ ^ 2 ( 3^2 - 1)(3^2 - 2) 

2 

3 . (1 - «1 - «2)(2 - - 3«2)(I - 3{i - 3«2) 

2 

^ _ 9 Si{ 2 ( 3 {i - 1 ) 

j,5 _f?2(3«2-l)(l-«l-«2) 

2 

. 9«l(l - «1 - e2)(2 - 3«l - 3«2) 

2 

9«l«2(3i2 - 1) 

N' = 

2 

h8 _ 9 ^ 2(1 -«1 -^2)(2 -3Ii - 342 ) 

2 

^9 _ 9Ii(3Ii-1)(1-4i-^2) 

2 

n10 = 27|i^2(1 -^1 -«2) 
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TABLE I. - CHARACTERISTICS OF THE STIFFNESS FINITE -ELEMENT MODELS 


DEVELOPED IN THE PRESENT STUDY 


Number of nodes dumber of internal : Total number of 

Element shape element degrees of freedom Approximation degrees of freedom, Designation 

per unknown 5r 


Quadrilateral 

4 

None 

Bilinear 

20 

SQ4 



1 

Bilinear plus bubble 
function 

25 

SQ5 


8 

None 

Quadratic 

40 

SQ8 



1 

Quadratic plus bubble 
function 

45 

SQ9 

Triangular 

6 

None 

Quadratic 

30 

ST6 


10 

1 

None 

Cubic 

50 

ST 10 


TABLE n.- THE NUMBER OF SYMBOLIC INTEGRATIONS NEEDED FOR TRIANGULAR AND QUADRILATERAL ELEMENTS 



Number of shape functions 
per element, r 

A*)*'' (l,j, 

k, { = 0 « r) 

(i, j = 0 - r; 

k = 1 - r; a = 1,2) 


O 

II 

r; 0 = 1,2) 

Element shape 

Number of numerically 
distinct integrals, 
(r+4)’./(4!r!)’ 

Number of representative 
integrals required 

Number of numerically 
distinct integrals 
r{r+l)(r+2) 

Number of representative 
integrals required 

Number of numerically 
distinct integrals, 
r(2r+l) 

Number of representative 
integrals required 

Triangular 

6 

210 

51 

336 

36 


78 

6 


10 

1001 

195 

1320 

121 


210 

13 

Quadrilateral 

4 

70 

17 

120 

11 


36 

3 


5 

126 

34 

210 

24 


55 

5 


8 

495 

84 

720 

54 

1 


136 

8 


9 

715 

130 

1 990 

1 83 

171 

1 1 


CO 





TABLE m. - CPU TIMES AND FIELD LENGTHS ON A CDC 6600 COMPUTER (USING THE RUN COMPILER 
UNDER SCOPE 3.0 OPERATING SYSTEM) FOR 14 SAMPLE PROBLEMS 


Element 

designation 

Element shape 

Required field 
len^hs* (octal) 

CPU time (in milliseconds) 
to compute A-, B-, and 
C -integrals 

Total CPU time 
(in milliseconds), 
plate elements** 

Total CPU time 
(in milliseconds), 
shell elements 

SQ4 

Parallelogram 

33 456 

8 

22 

36 

1 

Trapezoid 


12 

26 

40 1 

1 

Trapezium 


18 

32 

46 1 

SQ5 

Parallelogram 

34 276 

10 

30 

52 


Trapezoid 


16 

36 

58 

i 

Trapezium 


32 

52 

74 

ST6 

Triangle 

35 270 

12 

40 

96 

t 

00 

1 

Parallelogram 

40 167 

18 

62 

212 1 

i 

Trapezoid 

40 167 

44 

88 

238 

1 

1 

Trapezium 

42 662 

170 

214 

364 

SQ9 

Parallelogram 

42 077 

20 

76 

260 


Trapezoid 

42 077 

58 

114 

298 1 

1 


Trapezium 

44 572 

212 

268 

452 

ST 10 

Triangle 

44 257 

28 

98 

414 

♦Initial] 

ty 6000 octal words more are needed to accommodate the loader and load tables. 


**Includes bending-extensional coupling. 



I 


I ■■■■■II I II 


I II 


I III ■ III I III II II I 


I 1 1 mill II I 


TABLE IV.- AN ESTIMATE FOR THE NUMBER OF ARITHMETIC OPERATIONS 
REQUIRED TO GENERATE THE ELEMENT STIFFNESS MATRIX USING 
THE CONVENTIONAL NUMERICAL QUADRATURE PROCEDURE 

M 

[K]. y [C] [B,] 


L '=1 

Computation step 

Product of weighting coefficients and 
matrix [C], [C]g^g = Hg [C]g^g 

Product of [Cjj]g^g [Bflg^gj. 

Product of [Be]T^g([Cc] 

Summation over ju quadrature points 
Total 


Operations 

Additions 

Multiplications 

0 

(8)(8)/i 

(8)(8)(5r)/i 

(8)(8)(5r)p, 

(5r)(5r + l)(g)^ 

(5r)(5r + l)(g)^ 

(5r)(5r + 1) 

2 ^ 

0 

1 r/i(45r + 137) 

(l00r2 + 340r + 64) p 
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TABLE V. - NUMBER OF FLOATING-POINT OPERATIONS REQUIRED TO GENERATE 
THE ELEMENT CHARACTERISTIC ARRAYS USING THE SYMINSE PROGRAM 




Triangular 

Parallelogram 

(or 

Computation step 

Operations 

elements 

rectangular) elements 



ST6 

STIO 

SQ4 

SQ5 

SQ8 

SQ9 

Evaluations of integrals 

Additions 

515 

1 771 

221 

351 

1 029 

1 371 


Multiplications 

1306 

4 745 

524 

860 

2 669 

3 609 

Formation of linear combinations 

Additions 

6132 

27 830 

2030 

3045 

14 148 

17 685 

for stiffness kJj 

Multiplications 

8490 

38 510 

2804 

4204 

19 592 

24 488 

Formation of linear combinations 

Additions 

99 

265 

46 

61 

172 

199 

for arrays Kjj, P^, and 

Multiplications 

136 

366 

63 

87 

237 

281 

Sum of these three steps 

Additions 

6746 

29 866 

2297 

3457 

15 349 

19 255 


Multiplications 

9932 

43 621 

3391 

5151 

22 498 

28 378 


TABLE VI. - RATIOS OF CPU TIMES AND FLOATING-POINT ARITHMETIC 
OPERATION COUNTS FOR EVALUATING INTEGRALS AS FUNCTIONS 
OF FORMING THE CHARACTERISTIC ARRAYS FROM 
LINEAR COMBINATIONS OF THOSE INTEGRALS 


[All values refer to the SYMINSE program] 


Element designation 

Ratio of 
CPU times* 

Ratio of 

addition counts** 

Ratio of 

multiplication counts** 

Parallelogram SQ4 

0.286 

0.106 

0.183 

Parallelogram SQ5 

.238 

.113 

.200 

Triangle ST6 

.142 

1 

.00826 

.151 

Parallelogram SQ8 

.0928 

.0719 

.135 

Parallelogram SQ9 

.0833 

.0767 

.146 

Triangle STIO 

.0725 

.0630 

.122 


* Using data from table m, 
**Using data from table V. 
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